setwd("./bc_dataset")
library(multinichenetr)
library(Seurat)
library(Connectome)
library(ggplot2)
library(cowplot)
library(ComplexHeatmap)
library(SingleCellExperiment)
library(dplyr)
source('../DifferentialConnectome.wrapper.r', local=T)
dir.create(file.path(getwd(), "result_tables"), showWarnings = FALSE)

1 Load the full cohort

In this case, the original raw counts matrix will be utilized (discrete counts).

# Load data sce
sce_subset_bc<-readRDS("D:/Data/datasets/BreastCancer/data/sce_cohort1_updated.rds")

sce_subset_bc<-multinichenetr::alias_to_symbol_SCE(sce_subset_bc, "human") %>% 
               multinichenetr::makenames_SCE()
# Convert to Seurat Obj and assign discrete counts
seurat_subset_sce <- CreateSeuratObject(counts = SingleCellExperiment::counts(sce_subset_bc), meta.data = as.data.frame(colData(sce_subset_bc)))

# free RAM space
rm(sce_subset_bc)

Idents(seurat_subset_sce) <- seurat_subset_sce@meta.data$subType    
# Re-assign idents (prevents issues during connectomics analysis) 
seurat_subset_sce@meta.data$ident<-Idents(seurat_subset_sce)

table(Idents(seurat_subset_sce))
## 
## Endothelial_cell      Cancer_cell             CD8T             CD4T 
##             6256            53451            18271            19239 
##     Myeloid_cell               NK       Fibroblast      macrophages 
##             6502             2271            22571             6898 
##           T_cell           B_cell           CD4REG              gdT 
##             3093            11500             5798             2152 
##        Mast_cell              DCs        monocytes              pDC 
##             1073             2211              606               20
seurat_subset_sce@meta.data
# visualize the total number of cell types in the dataset
table(Idents(seurat_subset_sce))
## 
## Endothelial_cell      Cancer_cell             CD8T             CD4T 
##             6256            53451            18271            19239 
##     Myeloid_cell               NK       Fibroblast      macrophages 
##             6502             2271            22571             6898 
##           T_cell           B_cell           CD4REG              gdT 
##             3093            11500             5798             2152 
##        Mast_cell              DCs        monocytes              pDC 
##             1073             2211              606               20
table(filter(seurat_subset_sce@meta.data, expansion_timepoint=='OnNE')$ident)
## 
## Endothelial_cell      Cancer_cell             CD8T             CD4T 
##             2154            20755             3071             3999 
##     Myeloid_cell               NK       Fibroblast      macrophages 
##             2397              488            10754             1764 
##           T_cell           B_cell           CD4REG              gdT 
##              893             2719             1053              251 
##        Mast_cell              DCs        monocytes              pDC 
##              213              559              108                1

Last, we remove pDCs: in the downstream analysis of OnNE, the number of pDC is 1, which is not comaptible for allowing the differential connectome analysis.

not_pdc_ids<-row.names(filter(seurat_subset_sce@meta.data, ident!='pDC'))
seurat_subset_sce<-seurat_subset_sce[,not_pdc_ids]

 

2 Contrasts of interest

The analysis will focused on two aspects that can be highlighted by different contrasts:

  1. PDL1 communication in preE vs PreNE

  2. Chemotactic and (lymph)angiogenic interactions in E (on-pre) vs NE (on-pre)

 

2.1 PDL1 communication in preE vs PreNE

  • CD274 (PDL1) was reported to have higher expression in subpopulations of T cells in expanders patentients

  • CD274 interaction was predicted using CellPhoneDB in the contrasts between PreE-PreNE

  • crucial interplay IFNG-PDL1 signaling, which rapresents the disease hallmark

2.1.0.1 Differential connectomics analysis principles

Here, we approach the analysis via the application of the standard pipeline provided in the differential connectomics vignette.

In essence, this approach allows to compare two conditions in a single differential connectome (weighted, directed network that describes LR pairs as weighted edges and cell idents as vertices). In this specific case, the authors of Connectome provided the a specific scoring system, called perturbation score. In this case, the prioritization focuses on promoting the ligand receptor pairs for those senders and recievers cell identities that experience the highest absolute value of fold change variation in the test condition vs control, as described below:

 

In addition, the differential connectomics analysis also computes the the weight norm log fold change metrics. This metrics is the sum of the logarithm of the fold change expression between coondition and control for the ligand and the receptor in the LR pair.

The differential connectomics analysis output dataframe stores these two scores as “score” and “weight.norm.lfc” respectively.

 

 

2.1.1 Performing differential connectomics analysis of PreE vs PreNE

# Appliacation of the standard differential connectomics analysis, according to vignette preE-vs preNE
PreE_vs_PreNE<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "PreE", "PreNE")

# write results
openxlsx::write.xlsx(arrange(PreE_vs_PreNE, desc(weight.norm.lfc)), 
                 "result_tables/PreE_vs_PreNE_connectome_DB.xlsx", rowNames = FALSE)

2.1.1.1 Evaluation of the presence of inferred CD274-mediated in the output

filter(PreE_vs_PreNE, ligand %in%"CD274" | receptor%in%"CD274")

CD274 is missing among the results.

Possible reasons:

  • Not called DE during DEA for filtering of DifferentialConnectome

  • Not in Connectome prior knowledge LR DB (FANTOM5)

2.1.3 Differential connectomics analysis using NicheNet v2 DB

# Download LR DB of NicheNetV2
lr_network = readRDS(url("https://zenodo.org/record/7074291/files/lr_network_human_21122021.rds"))
dplyr::filter(lr_network, lr_network$from%in%"CD274" | lr_network$to%in%"CD274")

NicheNet V2 DB encodes multiple CD274-related LR pairs.

PreE_vs_PreNE_nichenet<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "PreE", "PreNE", db=as.data.frame(lr_network))


  
# Write xlsx based on desc weight norm lfc, not filtered
openxlsx::write.xlsx(arrange(PreE_vs_PreNE_nichenet, desc(weight.norm.lfc)),
                 "result_tables/PreE_vs_PreNE_connectome_via_nichenetLR_V2_DB.xlsx", rowNames = FALSE)
arrange(PreE_vs_PreNE_nichenet, desc(weight.norm.lfc))

2.1.4 Evaluation of CD274 pairs and conclusions

# Show CD274 LR reordered by weight.norm.lfc
PreE_vs_PreNE_nichenet%>% 
    arrange( desc( weight.norm.lfc))%>%
    mutate(rank_wnlfc= rank(desc(.$weight.norm.lfc)))%>%
    filter(ligand %in%"CD274" | receptor%in%"CD274")%>%
    select(6:18)

Conclusions: the analysis performed using NicheNet V2 prior model allows to capture and prioritize CD247 LR interactions. However, in the connectome paper, in the context of differential connectomics analysis the authors suggest:

(…)It is generally useful to limit analysis to those edges which have ligand and receptor expression in > 10% of their respective clusters in either the control or the test condition to avoid noisy measurements.; this thresholding functionality is built into CircosDiff and DifferentialScoringPlot.

This thresholding strategy is indeed applied in the mentioned functions as filtering criterias for these downstream visualization steps of the differential connectomics output (not yet applied because beyond of the scope of these markdowns). If we apply this to our final results, all the CD274 related LR interactions get filtered out.

filter(PreE_vs_PreNE_nichenet, pct.source.1 >0.1 &  pct.source.2  >0.1 & 
         pct.target.1  >0.1 & pct.target.2  >0.1 ) %>%
filter(ligand %in%"CD274" | receptor%in%"CD274")

2.2 Inference of Chemotactic and (lymph)angiogenic interactions

Connectome does not allow to accommodate an E(On vs pre) vs NE(on vs pre) contrast design directly. Therefore, we will compute and explore the outputs of the E(On vs Pre) and NE(On vs Pre) as workaround. In the following code we perform differential connectomics analysis via:

  • Differential connectomics analysis of OnE vs PreE and OnNE vs PreNE, standard Connectome vignette

  • Differential connectomics analysis of OnE vs PreE and OnNE vs PreNE via NicheNet V2 LR DB

# Appliacation of the standard differential connectomics analysis, according to vignette preE-vs preNE
OnE_vs_PreE<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "OnE", "PreE")
OnNE_vs_PreNE<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "OnNE", "PreNE")

# Appliacation of the standard differential connectomics analysis using NicheNet V2 LR DB
OnE_vs_PreE_nichenet_db<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "OnE", "PreE",  db=as.data.frame(lr_network))
OnNE_vs_PreNE_nichenet_db<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "OnNE", "PreNE",  db=as.data.frame(lr_network))


# write results, arranged by weight.norm.lfc
openxlsx::write.xlsx(arrange(OnE_vs_PreE, desc(weight.norm.lfc)), 
                 "result_tables/OnE_vs_PreE_connectome_DB.xlsx", rowNames = FALSE)

# write results, arranged by weight.norm.lfc
openxlsx::write.xlsx(arrange(OnNE_vs_PreNE, desc(weight.norm.lfc)),
                 "result_tables/OnNE_vs_PreNE_connectome_DB.xlsx", rowNames = FALSE)



openxlsx::write.xlsx(arrange(OnE_vs_PreE_nichenet_db, desc(OnE_vs_PreE_nichenet_db)), 
                 "result_tables/OnE_vs_PreE_connectome_via_nichenetLR_V2_DB.xlsx", rowNames = FALSE)

# write results, arranged by weight.norm.lfc
openxlsx::write.xlsx(arrange(OnNE_vs_PreNE_nichenet_db, desc(OnNE_vs_PreNE_nichenet_db)),
                 "result_tables/OnNE_vs_PreNE_connectome_via_nichenetLR_V2_DB.xlsx", rowNames = FALSE)
arrange(OnE_vs_PreE, desc(weight.norm.lfc))
arrange(OnNE_vs_PreNE, desc(weight.norm.lfc))
#NicheNet V2 DB Otuputs
arrange(OnE_vs_PreE_nichenet_db, desc(weight.norm.lfc))
arrange(OnNE_vs_PreNE_nichenet_db, desc(weight.norm.lfc))

2.2.1 Chemotactic interactions

From Bassels et al:

We next identified DEGs between T cells that expand versus T cells that do not expand pre-treatment (Fig. 4a, Extended Data Fig. 5a,b and Supplementary Dataset 5). Non-expanding T cells were more naive (LEF1, SELL, TCF7), while expanding T cells exhibited high effector function (IFNG), immune cell-homing signals (CXCL13, CCL3/4/5), cytotoxicity (GZMB, PRF1, NKG7), antigen presentation (CD74, HLA-DRB1/5, HLA-DQA2)

From MultiNicheNet paper:

Examples of increasing chemotactic interactions are CCL19-CCR7 between fibroblasts and B cells, CXCL12-CXCR4 between endothelial cells and B cells, and CXCL9-CXCR3 between macrophages/fibroblasts/DCs and (regulatory) CD4 T cells

chemotactic<-c('CXCL13', 'CCL3', 'CCL4', 'CCL5', 'CCL19', 'CXCL9' )

# check if genes are encoded in the connectome prior DB
chemotactic%in%connectome.genes
## [1] TRUE TRUE TRUE TRUE TRUE TRUE
OnE_vs_PreE%>%
  mutate(rank_wnlfc= rank(desc(.$weight.norm.lfc)))%>%
  filter(ligand%in%chemotactic)
OnE_vs_PreE_nichenet_db%>%arrange( desc( weight.norm.lfc))%>%
  mutate(rank_wnlfc= rank(desc(.$weight.norm.lfc)))%>%
  filter(ligand%in%chemotactic)

2.2.2 Lymphoangiogenic interactions

lymphangiogenesis<-c('VEGFC', 'FLT4', 'ITGA1', 'ITGA2', 'ITGA5')

# check if genes are encoded in the connectome prior DB
lymphangiogenesis%in%connectome.genes
## [1] TRUE TRUE TRUE TRUE TRUE
OnNE_vs_PreNE%>%
  mutate(rank_wnlfc= rank(desc(.$weight.norm.lfc)))%>%
  filter(ligand%in%lymphangiogenesis)
OnNE_vs_PreNE_nichenet_db%>%arrange( desc( weight.norm.lfc))%>%
  mutate(rank_wnlfc= rank(desc(.$weight.norm.lfc)))%>%
  filter(ligand%in%lymphangiogenesis)

Conclusions: no lymphoangiogenic interactions using NicheNet and mentioned in the paper were found using Connectome.

## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_Belgium.utf8  LC_CTYPE=English_Belgium.utf8   
## [3] LC_MONETARY=English_Belgium.utf8 LC_NUMERIC=C                    
## [5] LC_TIME=English_Belgium.utf8    
## 
## time zone: Europe/Brussels
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    grid      stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] gtools_3.9.4                plotrix_3.8-2              
##  [3] tibble_3.2.1                dplyr_1.1.3                
##  [5] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.2
##  [7] Biobase_2.60.0              GenomicRanges_1.52.0       
##  [9] GenomeInfoDb_1.36.4         IRanges_2.34.1             
## [11] S4Vectors_0.38.2            BiocGenerics_0.46.0        
## [13] MatrixGenerics_1.12.3       matrixStats_1.0.0          
## [15] ComplexHeatmap_2.16.0       cowplot_1.1.1              
## [17] ggplot2_3.4.3               Connectome_1.0.1           
## [19] SeuratObject_4.1.4          Seurat_4.4.0               
## [21] multinichenetr_1.0.3       
## 
## loaded via a namespace (and not attached):
##   [1] progress_1.2.2            nnet_7.3-19              
##   [3] locfdr_1.1-8              goftest_1.2-3            
##   [5] Biostrings_2.68.1         vctrs_0.6.3              
##   [7] spatstat.random_3.1-6     digest_0.6.33            
##   [9] png_0.1-8                 shape_1.4.6              
##  [11] proxy_0.4-27              ggrepel_0.9.3            
##  [13] deldir_1.0-9              parallelly_1.36.0        
##  [15] MASS_7.3-60               reshape2_1.4.4           
##  [17] httpuv_1.6.11             foreach_1.5.2            
##  [19] withr_2.5.1               xfun_0.40                
##  [21] ggpubr_0.6.0              ellipsis_0.3.2           
##  [23] survival_3.5-7            memoise_2.0.1            
##  [25] ggbeeswarm_0.7.2          emmeans_1.8.8            
##  [27] zoo_1.8-12                GlobalOptions_0.1.2      
##  [29] pbapply_1.7-2             Formula_1.2-5            
##  [31] prettyunits_1.2.0         KEGGREST_1.40.1          
##  [33] promises_1.2.1            httr_1.4.7               
##  [35] rstatix_0.7.2             globals_0.16.2           
##  [37] fitdistrplus_1.1-11       rstudioapi_0.15.0        
##  [39] miniUI_0.1.1.1            generics_0.1.3           
##  [41] base64enc_0.1-3           zlibbioc_1.46.0          
##  [43] ScaledMatrix_1.8.1        ggraph_2.1.0             
##  [45] polyclip_1.10-6           randomForest_4.7-1.1     
##  [47] GenomeInfoDbData_1.2.10   xtable_1.8-4             
##  [49] stringr_1.5.0             doParallel_1.0.17        
##  [51] evaluate_0.22             S4Arrays_1.0.6           
##  [53] hms_1.1.3                 irlba_2.3.5.1            
##  [55] colorspace_2.1-0          visNetwork_2.1.2         
##  [57] ROCR_1.0-11               reticulate_1.32.0.9002   
##  [59] spatstat.data_3.0-1       magrittr_2.0.3           
##  [61] lmtest_0.9-40             readr_2.1.4              
##  [63] later_1.3.1               viridis_0.6.4            
##  [65] lattice_0.21-9            spatstat.geom_3.2-5      
##  [67] future.apply_1.11.0       genefilter_1.82.1        
##  [69] XML_3.99-0.14             scattermore_1.2          
##  [71] scuttle_1.10.2            shadowtext_0.1.2         
##  [73] RcppAnnoy_0.0.21          class_7.3-22             
##  [75] Hmisc_5.1-1               pillar_1.9.0             
##  [77] nlme_3.1-163              iterators_1.0.14         
##  [79] caTools_1.18.2            compiler_4.3.1           
##  [81] beachmat_2.16.0           stringi_1.7.12           
##  [83] gower_1.0.1               tensor_1.5               
##  [85] minqa_1.2.6               lubridate_1.9.3          
##  [87] plyr_1.8.9                crayon_1.5.2             
##  [89] abind_1.4-5               scater_1.28.0            
##  [91] blme_1.0-5                locfit_1.5-9.8           
##  [93] sp_2.1-0                  bit_4.0.5                
##  [95] graphlayouts_1.0.1        UpSetR_1.4.0             
##  [97] codetools_0.2-19          recipes_1.0.8            
##  [99] BiocSingular_1.16.0       bslib_0.5.1              
## [101] e1071_1.7-13              GetoptLong_1.0.5         
## [103] plotly_4.10.2             remaCor_0.0.16           
## [105] mime_0.12                 splines_4.3.1            
## [107] circlize_0.4.16           Rcpp_1.0.11              
## [109] sparseMatrixStats_1.12.2  blob_1.2.4               
## [111] knitr_1.44                utf8_1.2.3               
## [113] clue_0.3-65               lme4_1.1-34              
## [115] listenv_0.9.0             checkmate_2.2.0          
## [117] DelayedMatrixStats_1.22.6 Rdpack_2.5               
## [119] openxlsx_4.2.5.2          ggsignif_0.6.4           
## [121] estimability_1.4.1        Matrix_1.6-1.1           
## [123] statmod_1.5.0             tzdb_0.4.0               
## [125] tweenr_2.0.2              pkgconfig_2.0.3          
## [127] tools_4.3.1               cachem_1.0.8             
## [129] RSQLite_2.3.1             RhpcBLASctl_0.23-42      
## [131] rbibutils_2.2.15          DBI_1.1.3                
## [133] viridisLite_0.4.2         numDeriv_2016.8-1.1      
## [135] fastmap_1.1.1             rmarkdown_2.25           
## [137] scales_1.2.1              ica_1.0-3                
## [139] nichenetr_2.0.4           broom_1.0.5              
## [141] sass_0.4.7                patchwork_1.1.3          
## [143] coda_0.19-4               carData_3.0-5            
## [145] RANN_2.6.1                rpart_4.1.21             
## [147] farver_2.1.1              aod_1.3.2                
## [149] tidygraph_1.2.3           mgcv_1.9-0               
## [151] yaml_2.3.7                DiagrammeR_1.0.10        
## [153] foreign_0.8-85            cli_3.6.1                
## [155] purrr_1.0.2               leiden_0.4.3             
## [157] lifecycle_1.0.3           caret_6.0-94             
## [159] uwot_0.1.16               glmmTMB_1.1.8            
## [161] mvtnorm_1.2-3             bluster_1.10.0           
## [163] lava_1.7.2.1              backports_1.4.1          
## [165] annotate_1.78.0           BiocParallel_1.34.2      
## [167] timechange_0.2.0          gtable_0.3.4             
## [169] rjson_0.2.21              ggridges_0.5.4           
## [171] progressr_0.14.0          parallel_4.3.1           
## [173] pROC_1.18.4               limma_3.56.2             
## [175] jsonlite_1.8.7            edgeR_3.42.4             
## [177] bitops_1.0-7              bit64_4.0.5              
## [179] Rtsne_0.16                spatstat.utils_3.0-3     
## [181] BiocNeighbors_1.18.0      zip_2.3.0                
## [183] muscat_1.14.0             jquerylib_0.1.4          
## [185] metapod_1.8.0             dqrng_0.3.1              
## [187] pbkrtest_0.5.2            timeDate_4022.108        
## [189] lazyeval_0.2.2            shiny_1.7.5              
## [191] htmltools_0.5.6           sctransform_0.4.0        
## [193] glue_1.6.2                factoextra_1.0.7         
## [195] XVector_0.40.0            RCurl_1.98-1.12          
## [197] scran_1.28.2              gridExtra_2.3            
## [199] EnvStats_2.8.1            boot_1.3-28.1            
## [201] igraph_1.5.1              variancePartition_1.30.2 
## [203] TMB_1.9.6                 R6_2.5.1                 
## [205] sva_3.48.0                tidyr_1.3.0              
## [207] DESeq2_1.40.2             gplots_3.1.3             
## [209] fdrtool_1.2.17            cluster_2.1.4            
## [211] ipred_0.9-14              nloptr_2.0.3             
## [213] DelayedArray_0.26.7       tidyselect_1.2.0         
## [215] vipor_0.4.5               htmlTable_2.4.1          
## [217] ggforce_0.4.1             car_3.1-2                
## [219] AnnotationDbi_1.62.2      future_1.33.0            
## [221] ModelMetrics_1.2.2.2      rsvd_1.0.5               
## [223] munsell_0.5.0             KernSmooth_2.23-22       
## [225] data.table_1.14.8         htmlwidgets_1.6.2        
## [227] RColorBrewer_1.1-3        rlang_1.1.1              
## [229] spatstat.sparse_3.0-2     spatstat.explore_3.2-3   
## [231] lmerTest_3.1-3            ggnewscale_0.4.9         
## [233] fansi_1.0.4               hardhat_1.3.0            
## [235] beeswarm_0.4.0            prodlim_2023.08.28